dtic  copy 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No.  0704-01-0188 

The  public  reporting  burden  for  this  collection  of  information  is  estimated  to  average  1  hour  per  response  including  the  time  for  reviewing  instructions,  searching  existing  data  sources, 
gathering  end  maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information  Send  comments  regarding  this  burden  estimete  or  eny  other  aspect  of  this  collection 
of  information,  including  suggestions  for  reducing  the  burden  to  Deportment  of  Defense,  Washington  Headquerters  Services  Directorate  for  Information  Operations  and  Reports 
(0704-0188)  1215  Jefferson  Devis  Highway,  Suite  1204.  Arlington  VA  22202-4302  Respondents  should  be  aware  that  notwithstanding  eny  other  provis  on  of  law  no  person  shall  be 
subject  to  any  penalty  for  failing  to  comply  with  a  collection  of  information  if  it  does  not  display  a  currently  valid  OMB  control  number 

PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADDRESS. 

1.  REPORT  DATE  (DD-MM-YYYY)  2.  REPORT  TYPE 

20-02-2009  REPRINT 

3.  DATES  COVERED  (From  -  To) 

4.  TITLE  AND  SUBTITLE 

Numerical  modeling  of  multidimensional  diffusion  in  the  radiation 
belts  using  layer  methods 


6.  AUTHORS 

Xin  Tao* 

Jay  M.  Albert 
Anthonv  A.  Chan* 


5a.  CONTRACT  NUMBER 


5b.  GRANT  NUMBER 


5C.  PROGRAM  ELEMENT  NUMBER 

62 1 0 1 OF 


5d  PROJECT  NUMBER 

1010 


5e.  TASK  NUMBER 

RS 


5f.  WORK  UNIT  NUMBER 
A 1 


7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

Air  Force  Research  Laboratory'  /RVBXR 

29  Randolph  Road 

Hanscom  AFB,  MA  01731-3010 


8.  PERFORMING  ORGANIZATION 
REPORT  NUMBER 

AFR1  ,-RV-H  A-TR-20 1 1-1016 


9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 


10.  SPONSOR/MONITOR’S  ACRONYM(S) 

AFRL/RVBXR 


11.  SPONSOR/MONITOR’S  REPORT 
NUMBER(S) 


12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  Public  Release;  distribution  unlimited 


13.  SUPPLEMENTARY  NOTES 

Reprinted  from  Journal  of  Geophysical  Research.  Vol.  I  14.  A022I5.  doi  10.1029/2008JA01 3826,  2009  ©2009.  American  Geophysical  Union 
*Department  of  Physics  and  Astronomy,  Rice  University,  Houston,  TX 


14.  ABSTRACT 

A  new  code  using  layer  methods  is  presented  to  solve  radiation  belt  diffusion  equations  and  is  used  to  explore  effects  of  cross 
diffusion  on  electron  fluxes.  Previous  results  indicate  that  numerical  problems  arise  when  solving  diffusion  equations  with  cross 
diffusion  when  using  simple  finite  difference  methods.  We  show  that  layer  methods,  which  are  based  on  stochastic  differential 
equations,  are  capable  of  solving  diffusion  equations  with  cross  diffusion  and  are  also  generalizable  to  three  dimensions.  We  run 
our  layer  code  using  two  chorus  wave  models  and  a  combined  magnetosonic  wave  and  hiss  wave  model  (MH  wave  model).  Both 
chorus  and  magnetosonic  waves  are  cqpable  of  accelerating  electrons  to  MeV  levels  in  about  a  day.  However,  for  the  chorus  wave 
models,  omitting  cross  diffusion  overestimates  fluxes  at  high  energies  and  small  pitch  angles,  while  for  the  MH  wave  model, 
ignoring  cross  diffusion  overestimates  fluxes  at  high  energies  and  large  pitch  angles.  These  results  show  that  cross  diffusion  is  not 
ignorable  and  should  be  included  when  calculating  radiation  belt  fluxes. 


15.  SUBJECT  TERMS 

Radiation  belt 


Cross  diffusion 


Electron  flux 


Chorus  wave 


Magnetosonic  wave 


16.  SECURITY  CLASSIFICATION  OF: 

17.  LIMITATION  OF 

18.  NUMBER 

19a.  NAME  OF  RESPONSIBLE  PERSON 

a.  REPORT 

b.  ABSTRACT 

C.  THIS  PAGE 

ABSTRACT 

OF 

PAGES 

10 

James  L  Metcalf 

UNCL 

UNCL 

UNCL 

UNL 

19B.  TELEPHONE  NUMBER  (Include  area  code) 

Standard  Form  298  (Rev.  8/98) 

Prescribed  by  ANSI  Std  Z39  18 


20110509004 


JOURNAL  OF  GEOPHYSICAL  RESEARCH,  VOL.  1 14,  A02215,  doi:  10. 1029/2008JA0 13826,  2009 


Click 

Her© 

for 

Fun 

Article 


Numerical  modeling  of  multidimensional  diffusion  in  the  radiation 
belts  using  layer  methods 

Xin  Tao,1  Jay  M.  Albert,2  and  Anthony  A.  Chan1 

Received  15  October  2008;  revised  26  November  2008;  accepted  17  December  2008,  published  20  February  2009. 

[i]  A  new  code  using  layer  methods  is  presented  to  solve  radiation  belt  diffusion  equations 
and  is  used  to  explore  effects  of  cross  diffusion  on  electron  fluxes.  Previous  results  indicate 
that  numerical  problems  arise  when  solving  diffusion  equations  with  cross  diffusion 
when  using  simple  finite  difference  methods.  We  show  that  layer  methods,  which  are  based 
on  stochastic  differential  equations,  are  capable  of  solving  diffusion  equations  with  cross 
diffusion  and  are  also  generalizable  to  three  dimensions.  We  run  our  layer  code  using  two 
chorus  wave  models  and  a  combined  magnetosonic  wave  and  hiss  wave  model  (MH  wave 
model).  Both  chorus  and  magnetosonic  waves  are  capable  of  accelerating  electrons  to 
MeV  levels  in  about  a  day.  However,  for  the  chorus  wave  models,  omitting  cross  diffusion 
overestimates  fluxes  at  high  energies  and  small  pitch  angles,  while  for  the  MH  wave  model, 
ignoring  cross  diffusion  overestimates  fluxes  at  high  energies  and  large  pitch  angles. 

These  results  show  that  cross  diffusion  is  not  ignorable  and  should  be  included  when 
calculating  radiation  belt  electron  fluxes. 

Citation:  Tao,  X.,  J.  M  Albert,  and  A.  A.  Chan  (2009),  Numerical  modeling  of  multidimensional  diffusion  in  lhe  radiation  belts 
using  layer  methods,  J  Geophys.  Res.,  114 ,  A02215,  doi:  10. 1029/2008JA0I3826. 


1.  Introduction 

[2]  The  Earth's  outer  radiation  belt  is  very  dynamic,  and 
electron  fluxes  can  vary  by  several  orders  of  magnitude  dur¬ 
ing  storm  times,  which  makes  it  very  hazardous  to  spacecraft 
and  astronauts  [e.g.,  Baker  et  al,  1986,  1994,  1997].  Quasi- 
linear  diffusion  theory  has  been  used  to  evaluate  dynamic 
changes  of  particle  fluxes  in  the  radiation  belts  [Albert,  2004; 
Albert  and  Youngs  2005;  Horne  and  Thorne ,  2003;  Horne 
et  al,  2003;  Li  et  al ,  2007;  Varotsou  et  al,  2005].  Albert 
[2004]  and  Albert  and  Young  [2005]  show  that  numerical 
problems  arise  when  solving  multidimensional  quasi-Iinear 
diffusion  equations  using  standard  finite  difference  methods, 
because  of  rapidly  v  arying  off-diagonal  terms.  Thus  different 
numerical  techniques  have  been  employed  to  solve  multidi¬ 
mensional  diffusion  equations,  c.g Albert  and  Young  [2005] 
and  Tao  et  al.  [2008]  (henceforth  denoted  as  paper  1). 

[3]  In  paper  1,  we  developed  a  stochastic  differential 
equation  (SDE)  code  to  solve  2-D  bounce-averaged  pitch 
angle  and  energy  diffusion  equations.  The  SDE  code  is  very 
efficient  when  solutions  on  a  small  number  of  points  are 
needed.  However,  if  solutions  are  needed  on  a  large  compu¬ 
tational  domain  for  long  times,  the  SDE  code  becomes  less 
efficient,  for  reasons  explained  in  paper  1 .  Milstein  [2002], 
Milstein  and  Tretyakov  [2002]  and  Milstein  and  Tretyakov 
[2001  ]  have  used  properties  of  numerical  integration  of  SDEs 
to  develop  so-called  layer  methods,  which  are  deterministic, 
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to  solve  parabolic  equations  successively  in  time.  In  this 
paper  wc  develop  a  code  using  layer  methods  and  show  that  it 
is  able  to  solve  2-D  radiation  belt  diffusion  equations  with 
cross  diffusion  and  it  is  generalizable  to  three  dimensions. 
Although  the  layer  code  does  not  have  the  high  paralleliza¬ 
tion  efficiency  compared  with  the  SDE  code  in  paper  1,  it  is 
more  efficient  when  solving  the  diffusion  equation  over  a  large 
computational  domain  for  long  times.  Also  our  layer  code 
can  handle  boundary  conditions  with  complicated  geometry 
rather  than  coordinate-equal-constant  type  boundaries  that 
are  typically  used  in  standard  finite  difference  codes. 

[4]  Using  this  layer  code,  we  then  explore  effects  of  ignor¬ 
ing  off-diagonal  terms  when  modeling  wave-particle  inter¬ 
actions  in  the  radiation  belts.  Previous  works  show  that  chorus 
waves  and  magnetosonic  waves  are  observed  during  storm 
times  and  are  possible  candidates  of  accelerating  electrons  to 
MeV  on  a  time  scale  of  days  [e.g.,  Horne  and  Thorne ,  1998; 
Horne  et  al,  2005,  2007;  Meredith  et  al,  2008].  Super- 
luminous  (RX,  LO,  LX  mode)  waves  have  also  been  shown 
to  be  possible  candidates  of  energizing  electrons  [Xiao  et  al , 
2006]  if  they  are  present  under  appropriate  conditions  [Xiao 
et  al ,  2007],  but  there  have  been  no  direct  observations  of 
these  waves  in  the  radiation  belts  so  far,  On  the  other  hand, 
interactions  with  EMIC  waves,  chorus  waves  and  hiss  waves 
have  been  invoked  as  important  loss  mechanisms  of  radiation 
belt  electrons  [e.g.,  Lyons  and  Thorne,  1973;  Summers  and 
Thorne,  2003;  Li  et  al,  2007].  In  this  work,  we  use  two  wave 
models:  the  chorus  wave  model  from  Li  et  al  [2007],  and  the 
combined  magnetosonic  wave  [Horne  et  al,  2007]  and  hiss 
wave  [Li  et  al,  2007]  (MH)  model.  In  paper  1,  we  showed 
that  ignoring  off-diagonal  terms  causes  errors  of  an  order  of 
magnitude  for  2  MeV  electrons  at  small  pitch  angles  using  the 
Horne  et  al  [2005]  chorus  wave  model.  Using  the  Li  et  al 
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[2007]  chorus  wave  model  and  the  MH  wave  model  is  helpful 
in  understanding  the  sensitivity  of  the  main  conclusions  of 
paper  1  to  different  wave  models. 

[5]  The  remainder  of  this  paper  is  organized  as  follows.  We 
introduce  the  layer  methods  by  using  a  simple  initial-value 
problem  in  seetion  2.  Details  of  our  2-D  layer  code  to  solve  a 
bounce-averaged  pitch  angle  and  energy  diffusion  equation 
are  given  in  seetion  3.  We  first  show  its  agreement  with  Albert 
and  Young  [2005]  results  in  section  3.1.  Then  we  solve  the 
diffusion  equation  with  diffusion  coefficients  calculated  using 
the  Li  et  al  [2007]  chorus  wave  model  (section  3.2)  and  the 
MH  wave  model  (section  3 .3 )  to  show  effects  of  ignoring  off- 
diagonal  terms  and  energization  of  electrons.  Our  results  are 
then  diseussed  and  summarized  in  section  4. 

2.  Milstein  Layer  Methods 

[6]  In  this  seetion,  we  will  use  an  initial-value  problem  to 
illustrate  the  layer  methods  shown  in  Milstein  [2002].  Bound¬ 
ary  conditions  can  be  implemented  in  a  similar  way  as 
described  by  paper  1,  or  by  Milstein  and  Tretyakov  [2001] 
and  Milstein  and  Tretyakov  [2002]. 

2.1.  One-Step  Representation  of  Solutions  Using 
the  SDE  Method 

[7]  Assume  that  we  want  to  solve  the  equation 


with  initial  condition//0,  x)  =  g(x).  First,  we  discretize  time  t 
equidistantly  to  /0,  •••,  /„,  /„+ 1,  and  assume  that  we  know 
solutions  of  all / at  time  t,n  which  means  that  now  //„,  x)  can 
be  considered  as  an  initial  condition  when  solving //„+!,  *). 
Then  using  the  SDE  method  described  in  paper  1,  we  have 

/(Om. (/(W))<  (2) 

where  S  is  the  expectation  value  and  x  is  given  by 

x  x+  b(tn+\,x)At  +  a(t„+\.x)AW.  (3) 

Here  At  =  tn+\  —  tn ,  a  is  related  to  a  by  <rcrr=  a  and  AW  = 
(A  W} ,  A  W2y . . .,  A  Wd)  is  one  increment  of  a  Wiener  random 
process  [ Gardiner ;  1985]. 

[8]  Numerically  AW  can  be  generated  from  a  vector  of 
standard  Gaussian  random  numbers  with  zero  mean  and  unit 
variance,  as  we  did  in  paper  1,  or  we  can  choose  the  prob¬ 
ability  distribution  of  the  components  to  be 


where  P  denotes  the  probability  [Milstein,  2002].  Substitut¬ 
ing  AW,  from  equation  (4)  into  equation  (3),  we  will  have  2d 
possible  Alf^s,  thus  2d  possible  x’s,  each  with  probability 
1/2^.  The  expectation  value  term  in  equation  (2)  can  then  be 
rewritten  as 

1  2" 

/(<.+,,*)*£  (/(/„,*))  =  =3  £/■('"  xj).  (5) 

y=i 


[9]  Note  that  in  reality,  we  usually  do  not  know / at  tn  and  at 
an  arbitrary  point  Xj.  This  is  why  we  trace  trajectories  baek  to 
the  initial  condition  in  paper  1.  However,  as  described  by 
Milstein  [2002],  we  can  use  interpolations  to  obtain //„,  x :j) 
from  already  known  fs  at  fixed  grid  points  to  make  a  con¬ 
vergent  algorithm  In  this  way,  we  obtain  solutions  succes¬ 
sively  from  time  layer  tn  to  /„+1,  hence  the  name  “layer 
methods”  [Milstein,  2002], 

2.2.  A  Simple  Layer  Method  Algorithm 

[10]  A  simple  interpolation  method  is  linear  interpolation. 
Take  a  1-D  diffusion  equation  for  example: 


d±-t& 

Ot  Ox  2  Ox 2 

(6) 

First,  discretizing  x  equidistantly  into  x0l  x2 

have 

XN,  WC 

/('">?)  =  v  /('"*')-  /('«.- 
X/+1  ~  X,  X,-!  —  X, 

+  0(  Ax2),  x,<x<xt,u 

(7) 

where  x,  and  jc/+1  are  fixed  grid  points,  and  Ax 
Then  a  simple  layer  method  algorithm  for  a  1  - 
problem  is 

=  X,+  \  -  xh 

-D  diffusion 

/frr+  i  Xj)  -^\f{tn,X j)  +/(4,X 2)]t  J  =  L2, 

...  ,M-  1, 

(8) 

Xi  =  Xj  4-  bAt  -+-  av^At 

(9) 

x2  =  Xj  +  bAt  —  <7\/A r, 

(10) 

=g(to.Xj), 

(11) 

with  //„,  x1<2)  calculated  using  equation  (7).  We  see  from 
equations  (7)  to  (11)  that  negative  values  of  f  cannot  arise 
from  this  procedure  A  proof  that  the  one-step  error  of  the 
above  algorithm  is  O(Ar)  is  given  in  Appendix  A.  Also  we 
show  the  connection  between  layer  methods  with  bilinear 
interpolation  and  conventional  finite  difference  methods  in 
Appendix  B  for  a  simple  2-D  diffusion  equation  without 
cross  diffusion. 

[11]  Dirichlet  and  Neumann  boundary  conditions  can  be 
implemented  in  a  similar  way  as  described  in  paper  1.  For 
example,  if  we  have  a  Dinehlet  boundary//,  x0)  =  g0(t ,  x0 ), 
we  replace  x  by  x0  if  x  <  x0 ;  if  we  have  a  Neumann  boundary 
d/Jdx  (/,  jcyv)  =  0,  then  we  replace  x  by  2x #  -  x  if  x  >  xN.  For 
more  general  ways  of  handling  boundary  conditions,  wc  refer 
readers  to  Milstein  and  Tretyakov  [2001, 2002] 

3.  Application 

[12]  In  this  seetion,  we  apply  the  layer  method  described 
in  section  2  to  the  bounced-averaged  equatorial  pitch  angle 
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(a0)  and  momentum  (p)  diffusion  equation  in  the  radiation 
belts 


±\/A /,  and  the  summation  in  equation  (18)  sums  over  four 
different  combinations  of  (A  IP),  AW2).  The  functions  b  and 
cr  are  the  same  as  in  paper  1 : 


2L 

dt 


11,/,  LJL+n 

Cp  dao  \  ‘ . p  dao  +  mp  dp, 

1  0  1  Of  „  0f\ 

*  Gdp° \Dtu,p p  da0  +  Dpp dp )  ’ 


(12) 


where  D«0n o,  D(t0p  and  Dpp  are  bounee-averaged  pitch  angle, 
mixed  and  momentum  diffusion  coefficients  [Albert,  2004]. 
Here  G  is  a  Jacobian  factor,  G  =  p27\ao)s\n(ao)cos(ao)y  and 
7\a0) «  1 .30  —  0.56sin(ao)  is  the  normalized  bounee  period. 
Initial  and  boundary  conditions  are  chosen  to  be  the  same  as 
those  of  Albert  and  Young  [2005]  and  paper  1 .  Thus  the  initial 
flux  is 


1  d  /gjW\  ,  1  d  fGD^p\ 
Gp  dao  \  p  )  G  Op  \  p  ) 


(21) 


<22> 

<?\\  =  x2DO0OJp  (23) 

*.2  =  0,  (24) 


j(t  =  0)  =  exp[— (£  —  0.2)/0. l][sin(a0)  -  sin(aar)],  (13) 

where  the  loss  eone  angle  q0l  =  5°  and  flux  j  is  related  to 
phase-space  density /by  j  =  J!p2.  Boundary  conditions  are 


^'loo=Oo/,  ^ ' 

(14) 

2L 

da  o 

=  0 

(15) 

f\t 

;.Emm=  0, 

(16) 

f  ~  0)  /Prntn 

(17) 

where  £min  =  0.2  MeV  and  Emax  =  5  MeV,  and  pmin  is  the 
momentum  corresponding  to  Emm  [Albert  and  Young ,  2005]. 

[13]  We  write  a  code  using  the  layer  method  to  solve 
the  diffusion  equation  (12).  Discretize  a0  and  y  =  log  E 
equidistantly  into  N„o  and  Ny  grid  cells,  and  thus  Aao  = 
(tt/2  -  a0L)/Na0  and  Ay  =  (ymax  -  YminV^y  The  equation  we 
use  to  solve /is 

/(r„+,.Qc  yj)=\YJt’" 

(18) 

where 

~  c*ot  +  banAt  +  (Ti,Alfi  +  cri2Alf2,  (19) 


Pj  -  Pj  4-  bpAi  +  a2\  AIT]  +  <t22 A W2,  (20) 

and  y,  is  then  obtained  from  pjt  If  a0/  <  Qol  or  ao,  >  rJ 2, 
we  replace  50/  by  a0L  or  it  -  50„  respectively.  If  /  <  pmin 
or  p  >  /W,  we  set  p  =  pmm  or  p  =  pmax,  respectively. 
In  equations  (18)  to  (20),  Atfri  and  A  W2  each  takes  value 


021  =  v  2  DaaP  /  JD  ooon  ■  (25) 


<ra  =  yJlDpp-a\v  (26) 

With  the  choice  of  bilinear  interpolation  to  obtain//*,  a0„>5/) 
in  equation  (18)  from  its  neighboring  grid  points  the  above 
algorithm  has  a  global  error  of  O(At)  when  Aa  =  cltAt,  Ay  = 
cvA/,  where  ca  and  cv  are  two  constants  [Milstein,  2002]. 
Milstein  [2002]  suggests  without  proof  that  it  may  be  possible 
to  use  eubic  interpolation  to  obtain  a  higher  convergence  rate 
We  conducted  numerical  experiments  using  eubie  interpola¬ 
tion  and  found  that  this  does  not  guarantee  positive  values 
Probably  a  more  sophisticated  method  would  give  both  high 
accuracy  and  positivity,  but  we  chose  to  use  bilinear  inter¬ 
polation  in  our  code  for  simplicity. 

3.1.  Comparison  With  Albert  and  Young  [2005]  Results 
[h]  To  show  that  layer  methods  can  be  used  to  solve  the 
diffusion  equation  (12),  wc  compare  results  calculated  using 
our  layer  code  with  results  of  Albert  and  Young  [2005]  using 
the  same  diffusion  coefficients.  The  diffusion  coefficients  are 
calculated  using  the  Horne  et  ai  [2005]  storm  time  ehorus 
wave  model  atL  =  4  5.  Wave  parameters  are  shown  in  Table  1 
for  a  eompanson  with  the  Li  et  al  [2007]  storm  time  chorus 
wave  model  (see  section  3.2). 

[15]  We  choose  Ar  =  4  x  10-4  day  to  give  a  relatively 
small  ehange  of  Aa0  and  Ay,  compared  with  the  computa¬ 
tional  domain.  We  plot  fluxes  of  0.5  MeV  and  2  MeV  at  /  =  2  d 
in  Figure  1  to  show  the  convergence  of  solutions  with  respect 
to  Nq0  and  Ny,  and  this  leads  to  our  choice  of  Na0  *  1 400  and 
Nv  =  1 500.  Comparisons  with  results  from  Albert  and  Young 
[2005]  are  shown  in  Figure  2  for  E  -  0,5  MeV  (top)  and 
E  =  2  MeV  (bottom)  electrons.  Considering  the  small  errors 
associated  with  eaeh  method,  we  conclude  that  the  two  sets  of 
results  agree  very  well  with  eaeh  other  and  our  layer  code  is 
capable  of  solving  the  bounce-averaged  pitch  angle  and 
momentum  diffusion  equation  (12)  with  cross  diffusion  In 
the  next  section,  we  apply  our  layer  code  to  the  diffusion 
equation  with  diffusion  coefficients  calculated  using  the  Li 
et  ai  [2007]  chorus  wave  model. 
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Table  1 .  Local  Time  Sector  Distribution,  Latitudinal  Distribution  |A|  of  the  Waves,  Equatorial  Ratio  of  Electron  Plasma  to  Gyro  Frequencies 
fpjfct*  and  Wave  Magnetic  Field  Amplitude  of  Wave  Models  Bw  From  Horne  ei  al.  [2005]  and  Li  et  al.  [2007] 


Home  et  al.  [20051 

Li  et  al. 

[20071 

2300-0600  MLT 

0600-1200  MLT 

1200-1500  MLT 

0000-0600  MLT 

0600  1200  MLT 

|A| 

0  to  15° 

15°  to  35° 

10°  to  35° 

0°  10  15° 

0°  to  35° 

f,Jf , 

B, 

3.43 

50  pT 

4.0 

100  pT 

6.72 

50  pT 

3.8 

50  pT 

4.6 

10o7$«o(ua  pT 

3.2.  Effects  of  Ignoring  Cross  Diffusion  in  Li  et  al  [2007] 
Chorus  Wave  Model 

[ 1 6]  Li  et  al  [2007]  used  a  new  chorus  wave  model  and 
calculated  changes  of  electron  fluxes  due  to  cyclotron  reso¬ 
nances  with  chorus  waves  by  solving  a  2-D  bounce-averaged 
pitch  angle  and  energy  diffusion  equation  using  an  implicit 
numerical  scheme.  However,  cross  diffusion  is  not  included 
in  their  calculation  [Li  et  al. ,  2007].  In  this  work,  wc  calculate 
diffusion  coefficients  including  cross  diffusion  using  the 
Li  et  al.  [2007]  main  phase  storm  time  chorus  wave  model. 
We  also  use  a  wave  normal  angle  distribution  from  Horne 
et  al  [2005],  in  contrast  to  Li  et  al  [2007],  who  use  a  parallel 
propagation  approximation.  The  resulting  diffusion  coeffi¬ 
cients  are  shown  in  Figure  3.  By  solving  the  diffusion  equa¬ 
tion  with  off-diagonal  terms  using  our  layer  code,  we  show  in 
this  section  effects  of  ignoring  off-diagonal  terms  on  electron 
fluxes  using  the  Li  et  al  [2007]  chorus  wave  model. 


Figure  1.  Fluxes  for  (top)  E  =  0.5  MeV  and  (bottom)  E  = 
2.0  MeV  at  t  -  2  d  with  different  choices  of  Nn0  and  N„ 
using  Horne  et  al  [2005]  chorus  wave  model.  In  the  top 
plot  the  three  lines  are  very  close  together,  in  contrast  to  the 
larger  separations  shown  in  the  bottom  plot. 


[n]  For  comparison  with  conclusions  in  paper  1 ,  the  bound¬ 
ary  and  initial  conditions  are  the  same  as  equations  ( 1 3 )-( 1 7), 
thus  they  are  different  from  those  in  Li  et  al  [2007].  Figure  4 
shows  color  plots  of  fluxes  calculated  using  the  diffusion 
coefficients  from  Li  et  al  [2007]  chorus  waves  at  t  =  0  1, 

1  and  2  days.  Both  results  with  (plots  on  left)  and  without 
(plots  on  right)  cross  diffusion  are  shown  for  comparison. 
We  see  from  Figure  4  that  at  high  energies  ignoring  cross 
diffusion  overestimates  fluxes  at  lower  pitch  angles  and 
creates  a  peak  in  fluxes  around  20°.  This  can  be  seen  more 
clearly  from  Figure  5,  which  shows  line  plots  of  fluxes 
calculated  with  and  without  cross  diffusion  at  t  =  0. 1 ,  1  and 

2  days  for  0.5  MeV  and  2  MeV  electrons.  At  t  =  0. 1  day,  the 
error  caused  by  ignoring  cross  diffusion  is  small  for  0  5  MeV 
electrons  at  all  pitch  angles  and  2  MeV  electrons  at  high  pitch 


Figure  2.  Comparisons  between  results  obtained  from  the 
layer  method  (solid  lines)  and  the  Albert  and  Young  [2005] 
method  (dashed  lines)  for  (top)  E  =  0.5  MeV  and  (bottom) 
£=20  MeV  at  t  =  0, 1  day  (blue  lines)  and  t  =  1 .0  day  (red 
lines)  Here  black  lines  show  the  initial  condition. 
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Figure  3.  Inverse  time  scales  in  units  of  s  1  from  diffu¬ 
sion  coefficients  calculated  using  Li  etal.  [2007]  chorus  wave 
model  with  the  wave  normal  angle  distribution  from  Horne 
et  al  [2005].  The  last  plot  shows  the  sign  of  the  cross  diffu¬ 
sion  coefficients. 


Figure  4.  Fluxes  calculated  by  the  layer  code  using  Li  et  ai 
[2007]  chorus  wave  model  at  t  =  0. 1 ,  1 ,  and  2  days.  The  plots 
on  the  left  show  fluxes  with  cross  diffusion,  and  the  plots  on 
the  right  show  fluxes  without  cross  diffusion. 


Figure  5.  Fluxes  for  (top)  E  —  0.5  MeV  and  (bottom)  E  = 
2.0  MeV  at  t  =  0.1  day  (blue  lines),  r  =  1  day  (green  lines), 
and  t  -  2.0  day  (red  lines)  with  and  without  off-diagonal 
diffusion  terms,  calculated  using  Li  et  al.  [2007]  chorus 
wave  model  Black  lines  show  the  initial  condition.  Dashed 
lines  are  results  without  off-diagonal  diffusion  coefficients, 
and  solid  lines  are  results  with  off-diagonal  terms. 


angles.  For  2  MeV  electrons  at  low  pitch  angles,  how  ever,  the 
error  is  about  a  factor  of  ~10.  At  /=  1  and  2  days,  at  0.5  MeV, 
ignoring  cross  diffusion  overestimates  fluxes  at  small  pitch 
angles  by  only  a  factor  of  2  ~  3.  However,  at  2  MeV,  ignor¬ 
ing  cross  diffusion  causes  an  error  of  about  two  orders  of 
magnitude  at  small  pitch  angles.  Thus,  similar  to  results  in 
paper  1,  ignoring  off-diagonal  terms  has  a  relatively  small 
effect  on  fluxes  for  lower-energy  electrons  at  higher  pitch 
angles,  but  it  introduces  larger  errors  for  larger  energy  elec¬ 
trons  at  lower  pitch  angles. 

[is]  To  understand  the  similarity  between  conclusions  ob¬ 
tained  here  and  in  paper  1 ,  we  now  discuss  features  of  the  Li 
et  al.  [2007]  and  Horne  et  ai  [2005]  wave  models,  whose 
parameters  are  listed  in  Table  1 .  First,  both  Li  etal.  [2007]  and 
Horne  et  al.  [2005]  wave  models  have  similar  latitudinal 
cutoffs  of  chorus  wave  power.  For  Li  et  al.  [2007]:  |  A|  <  35° 
on  the  daysidc,  and  A  <  1 5°  on  the  nightside,  for  Horne  et  al 
[2005]:  1 5°  <  |  A|  <  35°  in  the  prenoon  sector,  10°  <  |  A|  <  35° 
for  the  afternoon  sector,  and  |A|  <  15°  on  the  nightside.  Wc 
sec  larger  errors  at  smaller  pitch  angles  with  both  models. 
Second,  even  though  the  Li  et  ai  [2007]  chorus  wave  model 
has  a  dayside  wave  power  increasing  with  latitude,  which 
gives  a  more  abrupt  cutoff  at  the  maximum  latitude  than  the 
Horne  et  ai  [2005]  wave  model,  the  actual  values  of  the  wave 
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Figure  6.  Inverse  time  scales  in  units  of  s“ 1  from  diffusion 
coefficients  calculated  using  a  combined  magnetosonic  wave 
[Home  et  al. ,  2007]  and  hiss  wave  [Li  et  al .,  2007]  model. 
The  last  plot  shows  the  sign  of  the  cross  diffusion  coefficients. 


amplitude  are  not  very  different.  At  A  =  0°,  the  drift-averaged 
wave  amplitude  for  the  Li  etal.  [2007]  wave  model  is  6/24  x 
5.6  pT  +  6/24  x  50  pT  =  13.9  pT,  while  the  Home  et  al 
[2005]  model  gives  7/24  x  50  pT  =  14.6  pT.  At  A  =  35°,  the 
drift-averaged  wave  amplitude  for  the  Li  et  al  [2007]  model 
in  main  phase  is  6/24  x  141 .25  pT  +  6/24  x  50  pT  =  47.8  pT, 
while  the  Home  et  al.  [2005]  model  gives  6/24  x  100  pT  + 
3/24  x  50  pT  =  31.25  pT.  Thus  we  see  that  both  models 
assume  zero  amplitude  above  35°  latitude  and  have  compa¬ 
rable  wave  power  levels,  so  it  is  not  too  surprising  to  see 
similar  conclusions  on  ignoring  cross  diffusion  from  the  two 
wave  models. 

3.3.  Evolution  of  Electron  Fluxes  Using  a  Model 
of  Fast  Magnetosonic  Waves  and  Hiss 

[19]  Interactions  with  fast  magnetosonic  waves  have  been 
recently  suggested  by  Home  et  al  [2007]  to  be  a  possible 
important  acceleration  mechanism.  Because  these  interac¬ 
tions  typically  only  involve  the  Landau  resonance  ( n  =  0), 
coupling  of  diffusion  in  ao  and  p  is  expected  to  be  especially 
important  for  them.  For  the  wave  model  given  by  Horne  et  al 
[2007],  the  quasi-linear  diffusion  coefficients  of  the  magneto¬ 
sonic  waves  are  nonzero  only  over  a  limited  range  of  pitch 
angle  and  energy  [see  Albert ,  2008,  Figure  9].  Thus  we  com¬ 
bine  the  magnetosonic  wave  model  of  Home  et  al  [2007] 
outside  the  plasmasphere  with  the  main  phase  hiss  wave  model 
in  plumes  from  Li  et  al.  [2007]. 

[20]  The  MLT  averaged  difftision  coefficients  from  mag¬ 
netosonic  waves  (60%)  and  hiss  waves  (15%)  are  shown  in 
Figure  6.  A  similar  numerical  experiment  as  Figure  1  is  used 
to  determine  that  Ntt0  =  1400  and  Ny  =  1 500  is  necessary  to 
obtain  accurate  solutions.  The  resulting  evolution  of  electron 


fluxes  are  plotted  in  Figure  7  at  /  =  0. 1 ,  1  and  2  day.  We  see 
from  the  plots  on  the  left  (results  with  D^p)  that  the 
magnetosonic  waves  can  accelerate  electrons  to  MeVon  time 
scales  of  a  day,  and  the  fluxes  show  a  peak  around  55°, 
producing  a  butterfly  distribution,  at  high  energies.  Compar¬ 
ing  the  plots  on  the  left  with  the  plots  on  the  nght  (results 
without  DaQp\  we  see  that  ignoring  cross  diffusion  over¬ 
estimates  fluxes  at  larger  energies  and  larger  pitch  angles 
(>55°),  which  is  different  from  the  effects  using  the  Home 
et  al.  [2005]  and  Li  et  al.  [2007]  chorus  wave  models.  The 
effects  of  ignoring  cross  diffusion  can  be  seen  more  clearly 
in  Figure  8,  which  shows  fluxes  versus  equatorial  pitch  angle 
at  /  =  0.1,  1  and  2  days  for  0.5  and  2  MeV.  We  see  that  for 
0.5  MeV  electrons,  effects  of  ignoring  cross  diffusion  arc 
small  at  all  pitch  angles.  However,  for  2  MeV  electrons, 
ignoring  cross  diffusion  overestimates  fluxes  at  large  pitch 
angles  (>55°)  by  a  factor  of  5  ~  10  at  all  three  times,  and 
by  a  factor  of  ^5  at  small  pitch  angles  at  /  =  1  and  2  days. 

4.  Summary  and  Discussion 

[21]  In  this  work,  wc  introduce  the  layer  method,  which  is 
based  on  the  SDE  method  of  paper  1,  to  solve  multidimen¬ 
sional  radiation  belt  diffusion  equations.  Compared  writh  the 
SDE  method,  the  layer  method  is  deterministic  and  more 
efficient  when  solutions  on  a  large  computational  domain  are 
needed  for  long  times.  Compared  with  finite  difference 
methods,  the  layer  methods  are  less  efficient,  but  generalize 
to  three  dimensions  easily  and  are  able  to  handle  complicated 
boundary  geometries.  We  apply  the  layer  method  to  a 
bounce-averaged  pitch  angle  and  energy  diffusion  equation 
and  obtain  excellent  agreement  with  a  previous  method  [Albert 
and  Young ,  2005]  using  the  Horne  et  al.  [2005]  chorus  wave 


With  D,  ,  without 


20  40  60  00  20  40  60  00  £) 

'*0  °o  (orbitrory  units) 

Figure  7.  Fluxes  calculated  by  the  layer  code  using  the 
MH  wave  model  at  t  =  0. 1 ,  1 ,  and  2  days.  The  plots  on  the  left 
show  fluxes  with  cross  diffusion,  and  the  plots  on  the  right 
show  fluxes  without  cross  diffusion. 
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radiation  belt  diffusion  with  important  cross  diffusion  includ¬ 
ed.  However,  to  accurately  model  radiation  belt  dynamics, 
more  accurate  wave  models  and  initial  and  boundary  con¬ 
ditions  should  be  used  For  example,  comparing  Figures  4 
and  7,  it  is  easy  to  see  that  errors  could  be  large  if  inaccurate 
wave  models  are  used.  Further  work  should  be  done  with 
improved  wave  models  as  they  become  available.  Also,  ef¬ 
fects  of  different  initial  conditions  and  boundary  conditions, 
such  as  relativistic  kappa-type  functions  [Xiao  et  al.y  2008], 
should  be  considered  to  improve  radiation  belt  modeling. 

Appendix  A:  One-Step  Error  of  the  Layer  Method 

[25]  To  calculate  the  one-step  error  of  the  layer  method, 
consider  a  1-D  initial  value  problem, 

%=b{,'x)%+'2°2{,'x)U'  S(*.x)-g(*).  (Al) 

With  the  layer  method,  we  discretize  t  cquidistantly  into  to, 
*2»  •  •  •  with  time  step  h ,  and  the  approximate  solution  of 
equation  (Al)  at  (tk+\y  x)  is  given  by 

/*+  iM  =  \f  tk,x  +  b(tk+[,x)  +  o{tk+\.x)\fh 

+  ^/[»*  x  +  b{tk+\,x)  -  a(tk+i,x)Vh  (A2) 


Figure  8.  Fluxes  for  (top)  E  =  0.5  MeV  and  (bottom)  E  = 
2.0  McV  at  t  ~  0.1  day  (blue  lines),  t  =  1  day  (green  lines), 
and  t  =  2.0  day  (red  lines)  with  and  without  off-diagonal 
diffusion  terms,  calculated  using  the  MH  wave  model  Black 
lines  show  the  initial  conditions.  Dashed  lines  are  results 
without  off-diagonal  diffusion  coefficients,  and  solid  lines 
are  results  with  off-diagonal  terms. 


model.  We  show  that  our  layer  code  is  able  to  solve  multi¬ 
dimensional  diffusion  equations  with  cross  terms. 

[22]  We  then  use  the  layer  code  to  evaluate  effects  of  ig¬ 
noring  cross  diffusion  using  the  Li  et  al  [2007]  chorus  wave 
model,  as  a  comparison  with  the  Horne  et  al  [2005]  chorus 
wave  model  used  in  paper  1.  The  main  conclusion  is  similar 
to  paper  1;  that  is,  ignoring  off-diagonal  terms  produces 
larger  errors  at  smaller  pitch  angles  and  higher  energies.  We 
show  in  section  3.2  that  this  similarity  might  be  due  to  the  fact 
that  both  wave  models  have  a  latitudinal  cutoff  at  35°  and 
comparable  wave  power  levels. 

[23]  In  section  3.3,  we  show  evolution  of  electron  fluxes 
using  a  combined  magnetosonic  wave  [Horne  et  al.,  2005] 
and  hiss  wave  model  [Li  et  al. ,  2007].  We  show  that,  despite 
pitch  angle  scattering  by  hiss  waves,  electrons  are  energized 
to  MeV  in  2  days  of  simulation.  Ignoring  cross  diffusion 
overestimates  fluxes  at  larger  pitch  angles  and  higher  ener¬ 
gies,  in  contrast  to  the  effects  of  ignoring  cross  diffusion 
using  the  Horne  et  al.  [2005]  and  Li  et  al.  [2007]  chorus  wave 
models.  Overall,  we  conclude  that  cross  diffusion  terms  arc 
important  and  should  be  included  when  modeling  diffusion 
of  electrons  in  the  outer  radiation  belt. 

[24]  With  the  layer  method  and  other  methods  [ Tao  et  al , 
2008;  Albert  and  Young ,  2005],  it  is  now  possible  to  simulate 


[26]  To  show  the  error  offk+\(x),  we  expand  f{tk+\  -  h , 
x  +  b  ±  (7 yfh')  using  Taylor  expansion  and  the  resulting 
equation  is 


f(tk+)  -h.x  +  b  =  =/{/*+!,*)- A 


(W±^)|  +  i(^±26^)g 

(A3) 


Inserting  equation  (A3)  into  (A2)  yields 

-*(!-*!-  5 °(*!)  (A4) 


Using  equation  (Al),  equation  (A4)  becomes 

/hiW  =/(t|j)rO(/i2),  (A5) 

which  shows  that  the  approximation  fi^\{x)  differs  from  the 
true  solution  J[t^+i,  *)  by  a  term  proportional  to  h 2. 

[27]  In  practice,  we  use  interpolations  to  obtain  J[tk,  x  + 
b  ±  <j\/h)  in  equation  (A2)  from  fixed  grid  points,  denoted 
as /fo>  x  +  b  ±  cr yfh).  For  example,  using  linear  interpolation, 
we  have 

/^.i+6±av/ij  0{Xx2).  (A6) 
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Thus  with  Ax  =  cxh ,  where  cx  is  a  constant,  the  error  from 
interpolation  is  also  0{h2)  and  equation  (A5)  is  not  changed. 
Overall  the  one-step  error  of  the  above  layer  method  is  Q(h2). 
Milstein  [2002]  shows  that  the  global  error  (the  error  accu¬ 
mulated  from  time  to  to  f)  of  the  layer  method  is  0{h). 


Appendix  B:  Relationship  Between  Finite 
Difference  Methods  and  Layer  Methods 

[28]  Consider  a  2-D  diffusion  equation  with  constant  dif¬ 
fusion  coefficients  and  no  cross  diffusion, 


(Bl) 


According  to  the  layer  method,  the  updated  value  of/(x,  y) 
after  a  time  step  Ar  is  just  the  average  of / at  the  four  points 
fix  ±Lx,y±  Ly,  f)>  where  Lx  =  \/2DxxAt  and  Ly  =  y/2DyyAt 
according  to  the  presenption  in  section  2.2.  Use  bilinear 
interpolation  on  a  regular  grid  with  spacing  (Ax,  Ay),  and 
take  A t  small  enough  that  Lx  <  Ax  and  Lv  <  Ay  so  that  the 
points  (x  ±  Lxy  y  ±  Ly)  lie  within  the  neighboring  grid  cells. 

[29]  Then  the  value  of f  at  grid  point  ij  and  time 
is  given  by 


4/T  1  —  Vv/J-  +2(1  rx)ryftJ  +  \  +  rxr\'fi+\j+\ 

+  2rx (l  -ry)f?_Xj  +  4(1  -r,)(l  -  ry)f”  +  2rx(\  -  />)/£„ 
+  rxr/?_ ly_,  4-  2(1  -  rx)rJJ_l  +  (B2) 

where  rx  -  LJ  Ax  and  ry  =  LJAy.  Note  that  is 
guaranteed  to  be  positive  since  all  the  coefficients  of fnt±\j±\ 
are  nonnegativc. 

[30]  On  the  other  hand,  the  conventional  explicit  finite 
difference  scheme  for  equation  (Bl)  can  be  written  as 


4/T1  =  2 C^+1  +  2 ejp_v  +  4(1  -  c,  -  cv)f”  +  2c +  2 

(B3) 

for  time  step  Ar,  where  cx  =  2D„Af/(Ax)2  and  c,  = 
2DyyAt/(Ay)2. 

[31]  If  rx  =  cx  and  ry  —  the  two  schemes  come  into  close 
agreement.  Then  the  difference  between  the  two  expressions 
for  4 flf !  can  be  recognized  as  the  finite  difference  expression 
forrJt/y(Ax)2(Ay)2(5^/7^2x^2y).  The  corresponding  difference 
in  d fidt  is 

D"D"a'd§?-  <B4> 

This  extraneous  term  vanishes  as  At  —*  0. 

[32]  The  conditions  rx  =  cx  and  rv  =  cv  are  equivalent  to 
2DxxAt/(Ax)2  =  1  and  2DyvAt/(Ay)2  =  1  ,  respectively.  Then 
the  combination  2DxxAt!{  Ax)2  +  2DyvAf/(  Ay)2  has  the  value  2, 
while  the  CFL  stability  criterion  for  the  explicit  scheme  (for 
the  original  diffusion  equation)  requires  it  to  be  less  than  one. 
Thus  for  the  simple  case  discussed  above,  the  layer  method 
may  be  interpreted  as  an  explicit  scheme  run  at  an  unstably 
large  time  step,  but  stabilized  by  the  small  extra  term  of 
O(At),  which  is  the  same  order  as  the  error  of  both  methods. 


This  is  analogous  to  grid  diffusivity  terms  added  to  finite  dif¬ 
ference  schemes  in  the  Lax  method 
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